

#computes a^b mod n
def modexp(a, b, n):
	if( b==0 ):
		return 1;
	elif( b%2==0):
		return (modexp(a,b/2,n)**2)%n;
	else:
		return (modexp(a,b-1,n)*a)%n;
		


#deterministic for any n<341,550,071,728,321
def millerrabinPrimeTest(n):
	alist=[2,3,5,7,11,13,17];
	
	if( n<=17 and alist.count(n)>0):
		return True;	

	s=0;
	d=n-1;
	while( d%2 == 0):
		s+=1;
		d/=2;

	for a in alist:
		x=modexp(a,d,n);
		if( x==1 or x==n-1):
			continue;

		completed=True;
		for r in range(1, s):
			x=modexp(x,2,n);
			if x==1:
				return False;
			if x==n-1:
				completed=False;
				break;
		if completed:
			return False;

	return True;
	

is_prime=millerrabinPrimeTest;	

